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Abstract. - We observe a crossover from strong to weak chaos in the spatiotemporal evolution 
of multiple site excitations within disordered chains with cubic nonlinearity. Recent studies have 
shown that Anderson localization is destroyed, and the wave packet spreading is characterized by 
an asymptotic divergence of the second moment m2 in time (as t^'"^), due to weak chaos. In the 
present paper, we observe the existence of a qualitatively new dynamical regime of strong chaos, 
in which the second moment spreads even faster (as t^"), with a crossover to the asymptotic 
law of weak chaos at larger times. We analyze the pecularities of these spreading regimes and 
perform extensive numerical simulations over large times with ensemble averaging. A technique 
of local derivatives on logarithmic scales is developed in order to quantitatively visualize the slow 
crossover processes. 
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Wave propagation within random potentials is an in- 
terdisciplinary research field applicable to many diverse 
systems; regardless of their classical or quantum nature, 
the overall wave behavior provides common ground for 
understanding transport properties. One such property - 
theoretically predicted by Anderson |1 and since labeled 
"Anderson localization" (AL) - is a halt of wave propaga- 
tion due to exponentially localized normal modes (NMs) 
of the random potential. The significance of AL has been 
evidenced in the past decades by a bevy of experimental 
observation; including optics [2] , acoustics fs] , microwaves 
[4], and matter waves [s]. 

In many experimental situations, AL can be strongly 
altered by the appearance of nonlinearity within the po- 
tential. These nonlinearities can be induced via the non- 
linear Kerr effect in disordered photonic lattices |6l, or 
atomic Bose-Einstein condensate interactions (controlled 
by Feshbach resonances) in optical lattices [7\. 

The question of the interplay between disorder and non- 
linearity - how the two complement, frustrate, or reinforce 
each other - is thus of strong importance. The theoret- 
ical study of AL in random nonlinear lattices has been 
advanced using several approaches including the studies 
of transmission [s] and stationary solutions pi. Recent 
research in dynamics of wave spreading within nonlinear 
disordered media focuses on the spatiotemporal evolution 



of wave packets, debating the asymptotic spreading law 
against an eventual blockage 10 - 15 . Previous numerical 



studies show that the second moment of a wave packe t 
starting from a single site excitation grows as t^^^ ImIIis]. 
However when starting from a distributed single normal 
mode state, faster growth is reported, though not quanti- 
tatively assessed 14-16. Theoretical expectations range 



from spreading without limits to a slowdown and restora- 
tion of AL. Different spreading characteristics are also 



claimed to be i^/^ 
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and a two regime 
This letter aims 



case with i^/^ and asymptotic r^ 
to clarify some of these controversies. 

We first show that using estimates for average spacings 
and nonlinear frequency shifts, three different evolution 
regimes can be identified. In particular, in contrast to pre- 
vious results we obtain a new regime of strong chaos which 
is accessible by initial multiple site excitations, but not by 
single site excitations. Contrary to previous studies, we 
perform extensive ensemble averaging over 1000 disorder 
realizations. As a result, smooth functional dependencies 
of wave packet characteristics on time are obtained. Us- 
ing a technique of smoothing and local differentiation on 
logarithmic scales, we are able to observe the fast spread- 
ing regime of strong chaos, proceeded by the predicted 
crossover into the asymptotic regime of weak chaos, as ar- 
gued in ref. 17 1. No further saturation and slowing down 
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of the asymptotic spreading process is observed on the 
longest time scales of observation. 

We study two different Hamiltonian models. The first 
is the one-dimensional disordered nonlinear Schrodinger 
(DNLS) equation 
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in which e; is the onsite energy chosen uniformly 
from a [—W/2, W/2] random distribution. By 5t^i = 
d'HD/d{i^pi), the equations of motion are generated: 



iipi = eiipi + Pli'il'^i'i - ipi+i - ip, 
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The second model considered is oscillators on a quartic 
Klein-Gordon (KG) lattice, given as 
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where ui and pi are respectively the generalized coordi- 
nate/momentum on the site / with an energy of £i, and 
ii are the disordered potential strengths chosen uniformly 
in [1/2,3/2]. Likewise, dfui — —dT-Lx/dui generates the 
equations of motion 



"; 



-em 



+ v^i^i+i +UI-1 - 2ui) 
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The two models have only onsite cubic nonlinearity, but 
the specific methods discussed here can also be generally 
applied to other nonlinear models with other powers of 
nonlinearity or long range nonlinear dependencies. 

Similar to /3 in the DNLS case, the nonlinear control 
parameter for the KG model is the total energy, E = 
J2i ^i — 0- Both models conserve the total energy; the 
DNLS model also conserves the total norm S = X^JV'iP- 
As a practical note between the two models: the KG model 
is much more computationally friendly, with numerical in- 
tegration speeds of two orders of magnitude faster. For 
small amplitudes an approximate mapping, /35 ~ 3WE, 
from the KG model to the DNLS model exists 18 . Fur- 



ther analytics will, in general, be discussed in terms of the 
DNLS model, since it now straightforward to adapt sim- 
ilar results for the KG model using the aforementioned 
mapping. 

By neglecting the nonlinear terms, the DNLS model (IT]) 
reduces to the linear eigenvalue problem XAi — eiAi — 
(A;+i -I- Ai^i). This leads to a set of NM amplitudes, 
A^^i, with NM frequencies of X^ € [-W/2 - 2, W/2 + 2] in 
a spectrum width of A = 4 -f M^. The coefficient choice of 
1/{2W) in the KG model (Isl) allows a linear reduction to 
the same eigenvalue problem as for the DNLS, but with 
the values of e/ = W{ei - 1) and A^ = WujI - W - 2; 
in this case, A = 1 + A/W is the width of the squared 
eigenfrequency spectrum, ujI G [1/2,3/2 -I- A/W]. 

The NM asymptotic spatial decay is given by A^^j ^ 
g^-i/HK) -vyjiere ^[\v) is the localization length. It is ap- 
proximated [19] in the limit of weak disorder {W <^ 1) as 



£.{K) < C(0) ~ 96W-'^. The NM participation number 
P:^ ~ ^/J2i^ti characterizes the NM spatial extent. An 
average measure of this extent is the localization volume 
V, which is on the order of 3.3^(0) for weak disorder and 
unity in the limit of strong disorder [20 . The average 
frequency spacing of NMs within a localization volume is 
then d w A/V. The two frequency scales d < A are thus 
expected to determine the packet evolution details in the 
presence of nonlinearity. 

Nonlinearity induces an interaction between NMs. The 
variables 0^ = ^j A^^iipi determine the complex time- 
dependent amplitudes of the NMs. Since all NMs are ex- 
ponentially localized in space, each of them is effectively 
coupled to a finite number of neighbor modes, i.e. the 
interaction range is finite. However, the strength of this 
coupling is proportional to a characteristic norm density, 
n = \(f>\'^. The frequency shift due to the nonlinearity is 
then d ^ jSn. 

We track the normalized forms of the NM norm densi- 
ties, z^ = |0,.|VE;.I'^mP' for DNLS. The KG counterpart 
is normalized energy density distributions. These den- 
sities are sorted on the center-of-norm coordinate X^, = 
"Y^i lAf ^, and two measures in NM space used: the partic- 



ipation number P — l/J2i^^i 
tity of strongest excited sites 



^ which queries the quan- 
and the second moment 
1712 = J2ij(^ ~ '^)^^i' (where D = ^^^z^) which probes 
distances between a distribution's tail and center. The 
ratio of the two measures C = P"^ /'m2 (the compactness 
index [15| ) quantifies the sparsity of a packet - thermal- 
ized distributions have C ~ 3, while C ^ 3 indicates either 
very sparse packets, or partial self-trapping. 

We consider compact wave packets at t = spanning a 
width L centered in the lattice, such that within L there is 
a constant initial norm density of riin and a random phase 
at each site (outside the volume L the norm density is 
zero) . In the KG case, this equates to exciting each site in 
the width L with the same energy density, £ = E/L, i.e. 
initial momenta of pi = zL\/2£ with randomly assigned 
signs. 

For (3 = and L < V, the packet will extend over V 
during the time t ^ 2T:/d; after that, AL stops further 
spreading, with a lowered norm density n(r) « {ninL)/V. 
For L > V , the norm density will not change appreciably 
up to r, so n(r) « rii„. For /3 > 0, the nonlinear frequency 
shift should be compared with the average spacing d. If 
f3n{T) < d, most of the NMs are weakly interacting with 
each other; hence, this regime is dubbed "weak chaos". 
Once 6 = d, the NM frequency renormalization begins to 
allow some NMs to interact resonantly, i. e. strongly, with 
each other. If /3n(T) > d, almost all NMs in the packet are 
resonantly interacting. This regime will be coined "strong 
chaos". If (5 > A, a substantial part of the wave packet 
will be self-trapped. This is due to the occurrence of a 
wave packet frequency shifting out of resonance with the 
finite-width linear spectrum. In fact, partial self-trapping 
will occur already for 6 > 2 since at least some sites in 
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the packet may be tuned out of resonance. For a single- 
site excitation L = 1 the strong chaos regime shrinks to 
zero width and one is left only with either weak chaos or 
self-trapping fll''12''14''15]. The key distinction of the 
multiple-site excitations is that they can spread in 
the strong chaos regime. 

Fig IT] summarizes the predicted regimes in a paramet- 
ric space for the case L = V, in which lines represent 
the regime boundaries 5 = d and 5 = 2. The lower 
boundary is analytically found, via d = A/(3.3^(0)) with 
^(0) = 9614^"^ being the weak disorder estimate. More so- 
phisticated numerical estimates of d 20 1 yield only slight 
corrections ioiW > 6. It should be noted that the regime 
boundaries in figjl] are NOT sharp, rather there is some 
transitional width between the regimes. The weaker the 
strength of disorder, the larger the window of strong chaos. 
Inversely, ior W > % the strong chaos window closes al- 
most completely. Ideally, we should utilize the smallest 
possible value of W . Computational limits restrict this, 
so we choose a reference oiW = A. It is important to note 
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Fig. 1: (Color online) Parametric space of disorder, W , vs. the 
frequency shift induced by nonlinearity, 8, for the DNLS model. 
The KG analog is obtained by the small amplitude mapping 
£ — !> 3VF5 (see inset of fig[3|. Three spreading regimes are 
shown for dynamics dictated by: (i) weak chaos (pale blue), 
(ii) strong chaos (green), and (iii) the onset of self-trapping 
(pale red). The three circles show the initial numerical values 
used in fig[2] 

that 5 will be reduced in time, since a spreading wave 
packet increases in size and drops its norm (energy) den- 
sity. This gives the following interpretation of fig{l] given 
an initial norm density, the packet is in one of the three 
regimes (for example, the three circles in figll]). A packet 
launched in the weak chaos regime stays in this regime, 
while one launched in the strong chaos regime spreads to 
the point that it eventually crosses over into the asymp- 
totic regime of weak chaos at later times. 

In order to observe the crossover, we use L = 21 (which 



is approximately equal to V ior W = 4) in system sizes 
of 1000 - 2000 sites. For DNLS, an initial norm den- 
sity of riin = 1 is used, so that initially 5 ^ (3. Non- 
linearities (£ for KG) are chosen within the three spread- 
ing regimes (see fig[l]), respectively j3 g {0.04,0.72,3.6} 
and £ £ {0.01,0.2,0.75}. Ecis.([2]4]) are time evolved us- 
ing SABA-class split-step symplectic integration schemes 
15 21 , with time-steps of dt ^ 10^^ — 10^^ up to a max- 
imum t ^ 10^ — 10^. Energy conservations are accurate to 
< 0.1% (for discussion of numerical accuracy in symplectic 
integrators, please see ref. 



15). 
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Fig. 2: (Color online) Upper row: Average log of second mo- 
ments (inset: average compactness index) vs. log time for the 
DNLS (KG) on the left (right), for W = 4,L = 21. Col- 
ors/letters correspond the three different regimes: (i) weak 
chaos - (b)lue, /3 = 0.04 {£ = 0.01), (ii) strong chaos - (g)reen, 
/3 = 0.72 (£ = 0.2), (iii) self-trapping - (r)ed, 13 = 3.6 (£ = 
0.75). The respective lighter surrounding areas show one stan- 
dard deviation error. Dashed lines are to guide the eye to 
~ t^'^, while dotted-dashed guides for ~ i^". Lower row: 
Finite difference derivatives for the smoothed m2 data respec- 
tively from above curves. 

Ensemble averages over disorder were calculated for 
1000 realizations and are shown in figJ2^ (upper row). In 
the regime of weak chaos we find a subdiffusive growth 
of TO2 at large times according to 7712 ^ t"', a < 1, with 
a compactness index C w 3. Note that the subdiffusive 
growth is difficult to see initially in fig. [2] for two rea- 
sons. Firstly, the logarithmic scaling hides any small ini- 
tial growth, and secondly, there is a characteristic time 
scale for the packet to spread from its initial preparation. 
In the regime of strong chaos we observe a faster subd- 
iffusive growth of 7712, with an additional slowing down 
at larger times, as expected from the predicted crossover. 
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The compactness index is also ^ « 3, as in the weak chaos 
regime. FinaUy, in the regime of partial self-trapping 7712 
grows, but the compactness index C, decreases in time sub- 
stantially. This indicates that a part of the wave packet is 
arrested, and another part is spreading. 

In order to quantify our findings, we smooth (log 7712) 
with a locally weighted regression algorithm [22], and 
then apply a central finite difference to calculate the local 
derivative 

d (log 7772) 



a(logi) = 



dlogt 



(5) 



The outcome is plotted in the lower row in fig(2] 

In the weak chaos regime the exponent a{t) increases up 
to 1/3 and stays at this value for later times. In the strong 
chaos regime a{t) first rises up to 1/2, keeps this value for 
one decade, and then drops down, as predicted. Finally, 
in the self-trapping regime we observe an even larger rise 
of a(^). Additionally, we also performed numerics for W G 
{1,2,6} with respective initial packetwidths of L = V^ G 
{361,91,11}. Results are qualitatively similar to those 
shown in fig. [2] and thus omitted for graphical clarity. 

Following the analysis in ref. fl7], the modes inside the 
packet interact in a nonintegrable way leading to chaotic 
dynamics. The norm (energy) diffusion is characterized by 
a diffusion rate D - l3'^n'^{V{lin)f , where ViPn) sa 1 - 
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is the probability of packet modes to resonate 
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It is this resonance probability that largely dictates 
whether the chaos is strong or weak. From 7772 ~ 1/77^ and 
the diffusion equation 7712 '~^ Dt, one obtains an equation 
1/77^ ^ /3(1 — e~^"/'')i^/^ that determines the subdiffusive 
spreading crossover from the regime of strong chaos to 
that of weak chaos 



77l2 



/3tV2, 



f3n/d > 1 (strong chaos) 



^-2/3^4/3^1/3^ /377/d < 1 (weak chaos) 



(6) 



The impact of the strong chaos regime is seen in the 
resonance probability; V ^ 1 ii (in is sufficiently larger 
than d. Such a situation can be generated for packets with 
large enough /3n (or energy density £ for KG) in which 
every mode in the packet resonates, and the condition for 
strong chaos yields faster spreading, 7^2 ~ t^^'^. These 
predictions for strong chaos are then observed &tt ^ 10'^ — 
10'' (KG: 10'' — 10^) in figl2l time averages in these regions 
over the green curves yield a « 0.49 ± 0.01 (KG: 0.51 ± 
0.02). 

With spreading continuing in the strong chaos regime, 
the norm density in the packet will decrease, and eventu- 
ally /377 < d. Then a dynamical crossover occurs to the 
slower weak chaos subdiffusive spreading. This crossover 
spans logarithmic time scales. Nevertheless, in the green 
curves of figJ2] clear decay in a to values below 1/2 is 
observed. Fits of the decay further suggest a « 1/3 at 
t ^ 10'° — 10". The challenge remains to directly ob- 
serve saturation at times accessible in computational ex- 
periments. 



The duration of a = 1/2 (and thus when the crossover 
occurs) is largely dependent on how deep in the strong 
chaos regime the state is initially. Since the boundaries 
between different regimes are NOT sharp, but rather have 
some characteristic width, ideally one should utilize the 
smallest possible value of W . This is shown in figJ3] for 
the KG model. For W € {1, 2}, a long plateau at a = 1/2 
is clearly observed. For W G {4,6}, the initial energy 
density approaches one of the boundary lines and likely 
crosses into a boundary window, in which a < 1/2. 




Fig. 3: (Color online) Spreading behavior in the strong chaos 
regime for the KG model, with an initial energy density of 
£ = 0.1. The four curves are for the disorder strengths of: 
W = 1- (r)ed, W = 2 - (g)reen, W = 4- (o)range, W = 6- 
(b)lue. Inset: the KG analog of the DNLS parametric space, 
figfl] The four points correspond to the disorder strengths used 
in the main portion of the figure. 

In the regime of self-trapping, a good portion of the 
excitation remains highly localized, while the remainder 
spreads (red curves in fig[2|. Therefore, P does not 
grow significantly, but the second moment does. Con- 
sequently, C drops and is a good indicator of the degree 
of self-trapping. The time evolution of ^ for excitations 
in different regimes is shown in the insets of figj2j In 
the regimes of weak and strong chaos, if self-trapping is 
avoided, the compactness index at largest computational 
times is C ~ 2.85 ± 0.79 (KG: 2.74 ± 0.83), as seen in the 
blue and green curves of figj2j This means that the wave 
packet spreads, but remains thermalized {( ~ 3). For 
the self-trapping regime (red curves), the compactness in- 
dex asymptotically decreases to very small values. Note in 
fig[2]at intermediate times, there is transient growth where 
a > 1/2; nonetheless it remains subdiffusive (a < 1). At 
larger times, this overshoot decreases. This is presumably 
due to some self-trapped states which interact strongly 
with the spreading part of the packet and release their 
norm (energy) into the thermal cloud at some time. These 
more complicated scenaria are not yet quantitatively un- 
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derstood, and certainly remain for future exploration. 

This crossover can be expected to show up in measure- 
ments of the heat conductivity k at finite norm (energy) 
densities. According to the heat equation dT{x,t)/dt — 
{K/c)d'^T/dx'^ where T is the temperature and c the spe- 
cific heat. Therefore the heat conductivity is proportional 
to the diffusion rate k = cD. For small norm (energy) 
densities, heat is proportional to the densities. Therefore 
we expect that for /3n > d, i.e. in the regime of strong 
chaos, K ^ T^ (here T is the temperature). For small 
enough temperatures one crosses over into the regime of 
weak chaos, and consequently we expect k ^ T^ . 

Let us summarize. In the presence of nonlinearity 
within one-dimensional disordered systems, Anderson lo- 
calization is destroyed. In this Letter, we use a tech- 
nique of ensemble averaging and local derivatives on log- 
arithmic scales. In contrast to previous results for sin- 
gle site excitations, we find that multiple site excitations 
can evolve either in the asymptotic regime of weak chaos, 
or in an intermediate regime of strong chaos (excluding 
self-trapping for strong nonlinearities). In the weak chaos 
regime the second moment ?Ti2 grows subdiffusively as t^'^ . 
In the strong chaos regime subdiffusion is faster, yielding 
7TJ2 ^ f^'^, with a subsequently slow (on logarithmic time 
scale) crossover to an asymptotic weak chaos law. 
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